OHI British Columbia | OHI Science | Citation policy
OHI British Columbia | OHI Science | Citation policy
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'Figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(raster)
library(data.table)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
### includes library(tidyverse); library(stringr);
### dir_M points to ohi directory on Mazu; dir_O points to home dir on Mazu
dir_git <- '~/github/spp_risk_dists'
### goal specific folders and info
dir_setup <- file.path(dir_git, 'data_setup')
dir_data <- file.path(dir_git, 'data')
dir_o_anx <- file.path(dir_O, 'git-annex/spp_risk_dists')
source(file.path(dir_git, 'data_setup/api_fxns.R'))Compare MPAs to biodiversity intactness:
Examine correlation between risk and variance? A cell with medium risk and high variance indicates balance of low-risk species with some at risk species. A cell with medium risk and low variance indicates a more systemically at-risk cell.
mean_rast <- raster(file.path(dir_git, 'output', 'mean_risk_raster_010deg.tif'))
var_rast <- raster(file.path(dir_git, 'output', 'var_risk_raster_010deg.tif'))
n_spp_rast <- raster(file.path(dir_git, 'output', 'n_spp_risk_raster_010deg.tif'))
mean_var_df <- data.frame(mean_risk = values(mean_rast),
var_risk = values(var_rast),
n_spp = values(n_spp_rast)) %>%
filter(!is.na(mean_risk)) %>%
filter(n_spp >= 5)
mean_var_sample <- sample_n(mean_var_df, 100000)
ggplot(mean_var_sample, aes(x = mean_risk, y = var_risk)) +
ggtheme_plot() +
geom_point(aes(color = n_spp), alpha = .1) +
scale_color_gradientn(colors = c('grey90', 'lightblue', 'blue1', 'blue3', 'blue4', 'purple4'),
values = c(0, .01, .05, .1, .25, 1))
coef_var_sample <- mean_var_df %>%
filter(mean_risk != 0) %>%
filter(n_spp > 5) %>%
sample_n(100000) %>%
mutate(coef_var = sqrt(var_risk) / mean_risk,
inv_cv = 1 / coef_var)
ggplot(coef_var_sample, aes(x = mean_risk, y = inv_cv)) +
ggtheme_plot() +
geom_point(aes(color = n_spp), alpha = .1) +
scale_color_gradientn(colors = c('grey90', 'lightblue', 'blue1', 'blue3', 'blue4', 'purple4'),
values = c(0, .01, .05, .1, .25, 1))
Inverse coefficient of variation, i.e. mean risk / standard deviation of risk, will tend to be higher for high means that correspond to low variance. Low variance occurs when all species within a region are similar in risk; this indicates a system-wide moderate risk. High variance occurs when species risks differ from the central tendency; this indicates a system with low-risk species and high risk species, most likely skewed toward low risk as long as mean is closer to 0 than 1. High ICV probably is a good candidate for MPAs.
MPA coverage will be compared to non-MPA areas within EEZs. The .csv codes MPAs as their numeric IUCN category (i.e. Ia = Ib = 1, II = 2, …, VI = 6); areas that are specifically designated as no-take (that are not already in categories Ia, Ib, and II) are coded as 7; areas that are designated in the WDPA database but do not fall into one of these categories are coded as -1.
mean_rast <- raster(file.path(dir_git, 'output', 'mean_risk_raster_010deg.tif'))
var_rast <- raster(file.path(dir_git, 'output', 'var_risk_raster_010deg.tif'))
n_spp_rast <- raster(file.path(dir_git, 'output', 'n_spp_risk_raster_010deg.tif'))
n_threat_rast <- raster(file.path(dir_git, 'output', 'n_threat_raster_010deg.tif'))
trend_rast <- raster(file.path(dir_git, 'output', 'trend_raster_010deg.tif'))
eez_rast <- raster(file.path(dir_git, 'spatial', 'eez_rast_010_wgs84.tif'))
mpa_df <- read_csv(file.path(dir_git, 'spatial', 'wdpa_mpa_area_010_wgs84.csv'))
# tmp_rast <- calc(eez_rast, fun = function(x) x %in% c(1:212, 214:255))
### excludes high seas and antarctica
risk_df <- data.frame(cell_id = 1:length(values(mean_rast)),
eez = values(eez_rast),
mean_risk = values(mean_rast),
var_risk = values(var_rast),
n_spp = values(n_spp_rast),
n_threat = values(n_threat_rast),
mean_trend = values(trend_rast)) %>%
mutate(icv = mean_risk / sqrt(var_risk)) %>%
filter(!is.na(mean_risk)) %>%
filter(eez %in% c(1:212, 214:254)) ### 255 is disputed areas - remove...
risk_mpa_df <- risk_df %>%
left_join(mpa_df, by = 'cell_id') %>%
mutate(mpa_pct = ifelse(is.na(mpa_pct), 0, mpa_pct),
wdpa_category = ifelse(is.na(wdpa_category), 0, wdpa_category)) %>%
left_join(read_csv(file.path(dir_git, 'spatial/rgn_names.csv')), by = c('eez' = 'rgn_id'))no_take_df <- risk_mpa_df %>%
mutate(mpa = wdpa_category %in% c(1, 2, 7) & mpa_pct > .5,
protection = 'Ia, Ib, II')
conserv_df <- risk_mpa_df %>%
mutate(mpa = wdpa_category %in% c(1:4, 7) & mpa_pct > .5,
protection = 'I-IV')
protect_df <- risk_mpa_df %>%
mutate(mpa = wdpa_category != 0 & mpa_pct > .5,
protection = 'Any')
ggplot(no_take_df %>% sample_n(10000) %>% arrange(mpa),
aes(x = n_spp, y = mean_risk, color = mpa)) +
ggtheme_plot() +
geom_point(alpha = .2) +
scale_x_log10(breaks = c(1, 10, 100, 1000))
ggplot(no_take_df %>% sample_n(10000) %>% arrange(mpa),
aes(x = n_spp, y = mean_risk, color = r1_label)) +
ggtheme_plot() +
geom_point(alpha = .2) +
scale_x_log10(breaks = c(1, 10, 100, 1000))
ggplot(no_take_df, aes(x = mpa, y = mean_risk)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean risk vs mpa: no take (I, II, no-take)')
ggplot(no_take_df, aes(x = mpa, y = icv)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Inverse coefficient of variation vs mpa: no take (I, II, no-take)')
ggplot(conserv_df, aes(x = mpa, y = mean_risk)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean risk vs mpa: conservation (I-IV)')
ggplot(protect_df, aes(x = mpa, y = mean_risk)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean risk vs mpa: any protected area')
ggplot(no_take_df, aes(x = mpa, y = mean_trend)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean trend vs mpa: no take (I, II, no-take)')
ggplot(conserv_df, aes(x = mpa, y = mean_trend)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean trend vs mpa: conservation (I-IV)')
ggplot(protect_df, aes(x = mpa, y = mean_trend)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean trend vs mpa: any protected area')
t.test(no_take_df$mean_risk[no_take_df$mpa],
no_take_df$mean_risk[!no_take_df$mpa],
alternative = 'greater')
x <- lm(mean_risk ~ mpa + factor(r1_label), data = no_take_df)
summary(x)
t.test(conserv_df$mean_risk[conserv_df$mpa],
conserv_df$mean_risk[!conserv_df$mpa],
alternative = 'greater')
t.test(protect_df$mean_risk[protect_df$mpa],
protect_df$mean_risk[!protect_df$mpa],
alternative = 'greater')Repeat the above analysis but with range-rarity-weighted maps.
mean_rr_rast <- raster(file.path(dir_git, 'output', 'mean_rr_risk_raster_010deg.tif'))
# var_rr_rast <- raster(file.path(dir_git, 'output', 'var_rr_risk_raster_010deg.tif'))
sr_rr_rast <- raster(file.path(dir_git, 'output', 'sr_rr_risk_raster_010deg.tif'))
sr_rr_threat_rast <- raster(file.path(dir_git, 'output', 'sr_rr_threat_raster_010deg.tif'))
trend_rr_rast <- raster(file.path(dir_git, 'output', 'trend_rr_raster_010deg.tif'))
eez_rast <- raster(file.path(dir_git, 'spatial', 'eez_rast_010_wgs84.tif'))
mpa_df <- read_csv(file.path(dir_git, 'spatial', 'wdpa_mpa_area_010_wgs84.csv'))
# tmp_rast <- calc(eez_rast, fun = function(x) x %in% c(1:212, 214:255))
### excludes high seas and antarctica
risk_rr_df <- data.frame(cell_id = 1:length(values(mean_rast)),
eez = values(eez_rast),
mean_rr_risk = values(mean_rr_rast),
# var_risk = values(var_rr_rast),
sr_rr_risk = values(sr_rr_rast),
sr_rr_threat = values(sr_rr_threat_rast),
mean_rr_trend = values(trend_rr_rast)) %>%
filter(!is.na(mean_rr_risk)) %>%
filter(eez %in% c(1:212, 214:254)) %>% ### 255 is disputed territories
left_join(read_csv(file.path(dir_git, 'spatial/rgn_names.csv')), by = c('eez' = 'rgn_id'))
risk_mpa_rr_df <- risk_rr_df %>%
left_join(mpa_df, by = 'cell_id') %>%
mutate(mpa_pct = ifelse(is.na(mpa_pct), 0, mpa_pct),
wdpa_category = ifelse(is.na(wdpa_category), 0, wdpa_category))no_take_rr_df <- risk_mpa_rr_df %>%
mutate(mpa = wdpa_category %in% c(1, 2, 7) & mpa_pct > .5,
protection = 'Ia, Ib, II')
conserv_rr_df <- risk_mpa_rr_df %>%
mutate(mpa = wdpa_category %in% c(1:4, 7) & mpa_pct > .5,
protection = 'I-IV')
protect_rr_df <- risk_mpa_rr_df %>%
mutate(mpa = wdpa_category != 0 & mpa_pct > .5,
protection = 'Any')
ggplot(no_take_rr_df %>% sample_n(100000) %>% arrange(mpa),
aes(x = sr_rr_risk, y = mean_rr_risk, color = mpa)) +
ggtheme_plot() +
geom_point(alpha = .2) +
scale_x_log10(breaks = c(1, 10, 100, 1000))
ggplot(no_take_rr_df, aes(x = mpa, y = mean_rr_risk)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean risk vs mpa: no take (I, II, no-take)')
ggplot(conserv_rr_df, aes(x = mpa, y = mean_rr_risk)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean risk vs mpa: conservation (I-IV)')
ggplot(protect_rr_df, aes(x = mpa, y = mean_rr_risk)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean risk vs mpa: any protected area')
ggplot(no_take_rr_df, aes(x = mpa, y = mean_rr_trend)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean trend vs mpa: no take (I, II, no-take)')
ggplot(conserv_rr_df, aes(x = mpa, y = mean_rr_trend)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean trend vs mpa: conservation (I-IV)')
ggplot(protect_rr_df, aes(x = mpa, y = mean_rr_trend)) +
ggtheme_plot() +
geom_boxplot(outlier.alpha = .1, outlier.color = 'grey60') +
facet_wrap( ~ r1_label) +
labs(title = 'Mean trend vs mpa: any protected area')
There seem to be two philosophies of MPA priority: place MPAs in impacted areas to reduce pressure and prevent further degradatation, or place MPAs to protect pristine places to maintain their pristine state. This suggests that we may see a bimodal distribution of mean risk within MPA cells relative to the distribution of mean risk in general.
ggplot(no_take_df, aes(x = mean_risk, fill = mpa, color = mpa)) +
ggtheme_plot() +
geom_density(alpha = .3, size = .25) +
facet_wrap( ~ r1_label, scales = 'free_y') +
labs(title = 'dist of mean_risk inside and outside MPAs - no take')
ggplot(no_take_df, aes(x = icv, fill = mpa, color = mpa)) +
ggtheme_plot() +
geom_density(alpha = .3, size = .25) +
facet_wrap( ~ r1_label, scales = 'free_y') +
labs(title = 'dist of ICV inside and outside MPAs - no take')
ggplot(conserv_df, aes(x = mean_risk, fill = mpa, color = mpa)) +
ggtheme_plot() +
geom_density(alpha = .3, size = .25) +
facet_wrap( ~ r1_label, scales = 'free_y') +
labs(title = 'dist of mean_risk inside and outside MPAs - I-IV')
ggplot(protect_df, aes(x = mean_risk, fill = mpa, color = mpa)) +
ggtheme_plot() +
geom_density(alpha = .3, size = .25) +
facet_wrap( ~ r1_label, scales = 'free_y') +
labs(title = 'dist of mean_risk inside and outside MPAs - any MPA')
ggplot(no_take_rr_df, aes(x = mean_rr_risk, fill = mpa, color = mpa)) +
ggtheme_plot() +
geom_density(alpha = .3, size = .25) +
facet_wrap( ~ r1_label, scales = 'free_y') +
labs(title = 'dist of rr-weighted mean_risk inside and outside MPAs - no take')
ggplot(conserv_rr_df, aes(x = mean_rr_risk, fill = mpa, color = mpa)) +
ggtheme_plot() +
geom_density(alpha = .3, size = .25) +
facet_wrap( ~ r1_label, scales = 'free_y') +
labs(title = 'dist of rr-weighted mean_risk inside and outside MPAs - I-IV')
ggplot(protect_rr_df, aes(x = mean_rr_risk, fill = mpa, color = mpa)) +
ggtheme_plot() +
geom_density(alpha = .3, size = .25) +
facet_wrap( ~ r1_label, scales = 'free_y') +
labs(title = 'dist of rr-weighted mean_risk inside and outside MPAs - any MPA')